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Abstract 

Models and methods that are able to accurately and efficiently predict the flows of low-speed rarefied gases are in high 
demand, due to the increasing ability to manufacture devices at micro and nano scales. One such model and method is 
a Fokker-Planck approximation to the Boltzmann equation, which can be solved numerically by a stochastic particle 
method. The stochastic nature of this method leads to noisy estimates of the thermodynamic quantities one wishes to 
sample when the signal is small in comparison to the thermal velocity of the gas. Recently, Gorji et al have proposed 
a method which is able to greatly reduce the variance of the estimators, by creating a correlated stochastic process 
which acts as a control variate for the noisy estimates. However, there are potential difficulties involved when the 
geometry of the problem is complex, as the method requires the density to be solved for independently. 

Importance sampling is a variance reduction technique that has already been shown to successfully reduce the 
noise in direct simulation Monte Carlo calculations. In this paper we propose an importance sampling method for 
the Fokker-Planck stochastic particle scheme. The method requires minimal change to the original algorithm, and 
dramatically reduces the variance of the estimates. We test the importance sampling scheme on a homogeneous 
relaxation, planar Couette flow and a lid-driven-cavity flow, and find that our method is able to greatly reduce the 
noise of estimated quantities. Significantly, we find that as the characteristic speed of the flow decreases, the variance 
of the noisy estimators becomes independent of the characteristic speed. 
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1. Introduction 

Recent technological advances have resulted in manufacturing processes that have made possible the production 
of mechanical devices that operate on the scale of micro and nanometers JT]. Such technologies include lab-on-a- 
chip devices, micro-heat exchangers, gas chromatographers and micro-jet actuators for control in aerospace. At such 
small scales, the Navier-Stokes-Fourier (NSF) equations are no longer able to accurately model gas flows, due to 
the length scales of macroscopic gradients approaching the length of the molecules mean free path, A. This results 
in the existence of a region known as the Knudsen layer near solid wall boundaries where the gas is prevented from 
relaxing to thermodynamic-equilibrium, invalidating the assumption that locally the gas is close to thermal equilibrium 
required for the NSF equations to be valid. 

The Boltzmann equation is a mesoscopic model that is considered to provide the most accurate description of 
rarehed gases beyond Newton’s laws. Before the advent of such small scale technologies, rarehed gas flows’ largest 
application area was that of supersonic atmospheric flows, where the Mach number of flow, Ma > 1. Currently, the 
prevalent method for numerically approximating the solution to the Boltzmann equation in such regimes is a stochastic 
particle method called direct simulation Monte Carlo (DSMC) |I2 IS. Due to the stochastic nature of the method. 
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DSMC becomes very inefficient for low-speed flows. Typically the Mach number, Ma « 1 for flows within micro 
and nano technologies, and for a given level of statistical error, the computational costs of DSMC scale as Ma a. 
This results in very long computation times for such calculations, and methods which are able to efficiently solve for 
low speed flows are highly desirable. 

Currently, there are two methods that are able to greatly reduce the variance of the desired thermodynamic outputs 
of DSMC calculations. The first, low-variance DMSC (LVDSMC), works by adapting the DSMC collision routine 
to calculate the evolution of the deviation /j — f - fn from a Maxwellian distribution 0. In low speed flows 
the deviation from equilibrium is small, allowing for a dramatic decrease in the variance of samples. An alternative 
method, variance reduced DSMC (VRDSMC), is able to work without significant changes to the DSMC algorithm 0. 
The method relies on importance sampling, which allows the algorithm to sample from an equilibrium distribution 
where the thermodynamic variables are known aproiri, to create estimators with smaller variance. 

More recently an alternative method to DSMC, where the Boltzmann collision operator is approximated by a 
Fokker-Planck operator, has been developed and shown to be more efficient than the basic DSMC algorithm 121 0. 
Like DSMC, it is solved stochastically using notional particles that represent a certain number of real particles in the 
gas to be simulated, and as such, the basic algorithm suffers from the same inhibitive scaling with the Mach number. 
Recently Gorji et al. El have proposed a method to reduce the variance of the Fokker-Planck solution algorithm 
that relies on creating a correlated equilibrium solution using the same set of random numbers that are used in the 
stochastic solution of the non-equilibrium process. The parallel correlated equilibrium process is in effect a control 
variate for the non-equilibrium process. 

This purpose of this paper is to develop an importance sampling variance reduction scheme for the Fokker-Planck 
method and demonstrate its effectiveness in simple test cases. The paper is organised in the following way: in section 
2 we introduce the Fokker-Planck model, and numerical stochastic particle scheme which we would like to create 
variance reduced estimators for. We then outline the general method that allows one to create variance reduced 
estimators by exploiting known information about how the macroscopic fields behave at equilibrium. In section 3 we 
describe the variance reduction scheme proposed by Goiji el al. Q, which creates a correlated equilibrium scheme. 
In section 4 we propose our importance sampling scheme, which we test in section 5 on a homogeneous relaxation, 
Couette flow and a lid-driven cavity flow. We then compare the importance sampling method against the results 
obtained by using a correlated equilibrium solution. 

2. Background 

2.1. The Fokker-Planck collision operator 

The Fokker Planck collision operator has appeared in several different contexts, originally derived for the distribu¬ 
tion function of a Brownian particle in a fluid iflOl . but can also be derived from an expansion of a linear Boltzmann 
equation, when considering the evolution of density function for a particle in a heat bath ini. It has been used to 
model electrons, dense liquids and more recently has received attention for its ability to model rarefied gas flows El- 
More recently it has been extended to describe flows of monatomic gas mixtures El, diatomic molecules 03 and 
has been coupled to DSMC El- The equation for the one-particle distribution function /(x, v, f), over a state-space 
comprised of the position x e velocity v e and time f 6 takes the form: 

+ V - V,/ = m{f) (1) 

:= ivvfcZ + RrVv/], (2) 

where t is a relaxation time, c = v - u is the local relative molecular velocity, u is the mean velocity: 

u(x, 0^ v/(x, V, t) dv, (3) 

T is the local temperature given by 

T(x,t) ^ ^ J c^f(x, V, f) dv, (4) 
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where R is the specific gas constant and p is the local density given by 


P(x, 0 



x,v, 0 dv. 


The collision operator has the property of conserving mass, momentum and energy. That is 


( 5 ) 


J^(/)^ufv = 0, (6) 

where if/ - v^] is the set of collisional invariants. The advantage of having a collision operator which can be writ¬ 

ten as a Fokker-Planck equation is that there exists an equivalent stochastic differential equation (SDE) representation 
for the dynamics of a random variable {X,, Vf) whose distribution / evolves according to Q; 


dX, = V,df (7) 

1 I '2.RT 

dV, = - (V, - U)df -H J — dW„ (8) 

where W, is a 3-dimensional Wiener process. An efficient scheme for evolving a collection of representative particles 
with positions and velocities j - \. .N, distributed according to the distribution f{x,\,t) in time was 

devised by Jenny et al. Q, and can be summarised as: 


Vi (f + Af) = y,.(f) - (1 - (y,(0 - Ui(t)) + - ^^24 (9) 

Xiit + At) = Xi(t) + Ui{t)At + T(y,(f) - (1 - (10) 

where i - 1,2,3 indexes the dimension, 

A ^ RT {l - (11) 

C = B^T(l-e^'/^)^ (13) 


and At is the time-step, r is the relaxation time and ^ are standard normal distributed random variables. The spatial 
domain is gridded into cells, and expectations of macroscopic quantities of interest are calculated during each time- 
step for each computational cell. The correct viscosity is obtained by choosing the relaxation time t = 2/i/p, where p 
is the viscosity and p is the pressure. 

2.2. Variance reduction for Monte Carlo sampling 

In this section we outline the general framework which allows one to reduce the variance of an estimator of a 
particular random variable, common to control variate and importance sampling schemes. Essentially, both methods 
exploit information about errors in estimates of known quantities. The general principal is as follows. Suppose we 
have a random variable X, and we wish to estimate E[X], let our estimate of E[X] be donated by X. Let T be a different 
random variable with known expectation E[T] with an estimator denoted by Y. Then for any cr 6 K we can use the 
identity 



E[X] = E[X-HaT] -q'E[T], 


(14) 
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to create a new unbiased estimator for E[X], 


XvR ^X + aY -aE[Y]. 


(15) 


The variance of this estimator is 


Var[lv«] = Var[X] + cp-V 2 iX{Y] + 2aCov{X, Y] (16) 

and if we minimise this over possible choices of a, then minimiser (^‘is given by 


* Coy[X,Y] 

a -7—, 

Var[F] 


(17) 


hence the variance for this choice of a is 


CoviX, 

Var[XyR] = Var[X]-. (18) 

Var[F] 

The only condition required for the variance of the estimator to be less than the variance of the original estimator 
is for Cov[X, T] > 0, and so X and Y being dependent is a necessary condition. This is all supposing that we already 
know, a* which presupposes that we already know Cov[X, T]. In general this is not something that is not known a 
priori, but can be estimated throughout the simulation. In the next sections we will see how in practice it is possible 
to exploit this. 

3. Parallel process variance reduction 

We now briefly describe the method proposed by Q. The objective of the method is to create a stochastic process 
Zf that is able to run in parallel to the original particle scheme, where crucially, the macroscopic fields are already 
known. If this is performed in a manner where the parallel process is correlated with the original stochastic process 
then the variance of the estimators can be reduced in the way described in the previous chapter. The coupling of the 
new stochastic process Z, to the original SDEs ([^ is achieved in the following way: 


dX, 

= V,df 

(19) 


1 llRT 


dV, 

= -(V,-u)df + J -dW, 

(20) 


T V T 


dZ, 

= Adf+ DdW, 

(21) 


where the coefficients A and D are chosen to keep the marginal distribution of (X,, Z,), which we denote /o(x, z, f), as 
a solution to a Fokker-Planck equation 


^ + V, ■ (U/o) = Vv ■ [a/o + ^Vv/o], 




dt * ' ' l 2 

which when supplemented by appropriate boundary conditions, is solved by a Maxwellian 


( 22 ) 


/m(x, z, 0 = ^ exp 

(InRTof'^ 


2RTo 


( 23 ) 
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The way to choose A and D to ensure that (23 1 is a solution of (22 1 is discussed in depth in i). The coupling of 
the processes Z, and V,, by the same Wiener process W, requires that when we use particles to generate numerical 
solutions, the same set of random numbers is used for both the equilibrium and non-equilibrium processes. This 
results in the correlation of estimations of expected quantities, allowing the kind of variance reduction outlined in the 
previous section. In this method, the parallel equilibrium process, with distribution /m, shares the same density p as 
the non-equilibrium process described which has / as its distribution function, and so the method cannot reduce the 
noise in the density calculation directly. Goiji et al propose that the density p is found from the continuity equation 
using conventional finite difference methods. The results they obtain for a homogeneous relaxation, Poisuille flow 
and lid-driven cavity flows show the method has the ability to substantially reduce noise. Because this method uses 
common random numbers to reduce the variance, we will refer to this method as a common random numbers (CRN) 
scheme. 


4. Importance sampling 

The method we propose is an importance sampling scheme. It differs from the importance sampling scheme 
utilised by VRDMSC 0, because the DSMC and Fokker-Planck method account for collisions in different ways. The 
principle that underpins the importance sampling remains the same however. Suppose we are interested in evaluating 
the expectation of g(v) where v is distributed according to the distribution function /, then given A independent 
samples {Vi,..., Va^) distributed according to / the following definition gives rise to the estimate; 


E/[.g(v)] 


N 




(24) 

(25) 


which we know from the Central Limit Theorem, has an error of order N We now define a weight function 


mv) := (26) 

/(v) 

which is a measure of how likely one is to see a particle with velocity v, relative to how likely one is to observe this 
particle if it was distributed to a reference density /ref. This definition is well defined if the distribution / is absolutely 
continuous with respect to /ef, meaning that frei{S ) = 0 whenever f{S ) = 0 for any subset S of the state-space. This 
definition can be viewed as a Radon Nikodym derivative. It can then be observed that the expectation of g(v) with 
respect to the reference distribution can be estimated using the original samples: 


E/,er[Mv)] = J fwf(y)g(y) dv 

(27) 


(28) 

= J f(y)W(v)g(v)dv 

(29) 

1 ^ 

^ - Xi ^(V')^(V'). 

/=1 

(30) 

This is significant as it allows one to sample from the reference distribution /ef, using the original set of samples 
from the distribution /. If the reference distribution is Maxwellian, /ef = /m where one knows the thermodynamic 
fields analytically, then one has the ability to create variance reduced estimators as described in section 2. In order to 
practically apply this, we need a method of evolving the weights and velocities {V', W] in time, where W‘ - W(V'). 
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For VRDSMC this is possible because it can be shown directly from the Boltzmann equation, that if two particles are 
chosen to collide with weights Wi and Wj then the post collision weights must be equal to ^(W + W-'). Because the 
Fokker-Planck dynamics have no explicit collisions, a different way to update the weights is needed. 

Importance weights can be initialised exactly, because the initial velocities of the particles are distributed according 
to a prescribed initial distribution /q. As time is evolved during the calculation, the distribution of velocities will 
change and hence so must the weights attached to each particle. VRDSMC is able to do this by creating collision 
rules that ensure that post collision velocities are still able to sample from the same reference distribution. However, 
these rules are not relevant for the Fokker-Planck particle dynamics as there are no explicit collisions. 

4.1. Weight update rule 

Instead, let us suppose that a given particle updates its velocity from V, ^ V,+i, where V? is distributed according 
to ft and V,+i is distributed according to f+i, and that we know Wt - VF(Vf). In order to update the weight exactly, 
one would need to know f+ilYt+i), however this distribution function is unknown. A simple method to estimate the 
updated weight is to use the zeroth order Taylor expansions of the joint distributions of V? and V,+i, resulting in the 
estimate: 


IV,+1 ^ IV,+1 


^(Vt+i|V,)/,,(¥,) 


/,rl(Vt+i|V,)/,(V,) 


^(Vt+i|V,) 

/,+i(Vt+i|V,) 


1V,(V,). 


(31) 

(32) 


This has approximation immediately has some desirable properties. Firstly, the error of the approximation decays 
with At. Also, it is possible to calculate this explicitly from the update rule V, ^ V,+i given by equation (j^. This 
conditional distribution will be a gaussian centred on V, plus the deterministic drift, with a temperature dependent 
variance. Further to this, it has the correct conditional expectation E[VV+i|lV,] = Wt when the distribution is stationary. 
However, on its own it is not a suitable choice as if such a rule is repeated the variance of this approximation diverges, 
which is a common problem for this type of particle weight importance sampling method Q Qh). This is a problem, 
because to reduce the variance of our estimators in a meaningful way, we require the weights to be close to unity. To 
avoid this problem we use the same kernel density estimator approach as used by the VRDSMC method |[6l. Kernel 
density estimation (KDE) is a method that allows one to obtain an estimate / of a density function / from samples 
distributed according to that density function in the following way: 


/(v)= ^2 k,(v-v'), (33) 

/=! 

where Kr is a kernel function that integrates over the state-space to 1, and r is a smoothing parameter that controls the 
width of the kernel function. We use the same spherical kernels as ISl: 


Kr(y - v') 


3/(47rr3) if||v-vi|<r 

0 otherwise 


(34) 


which returns the reciprocal of the volume of a sphere of radius r if v' lies within the sphere of radius r centred on v, 
and otherwise returns a zero. If we combine this with ([3^, the update rule that is obtained is 


lV,+i(V') 


SjLi ^r(Vi - V-')lV,+l(V^) 

zL^xVi-v^) 


\Srm\ 




VJeSrW) 


(35) 

(36) 
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where 5r(V') = {V^ : HV-' - V‘|| < r} is the set of samples whose members lie within the sphere of radius of r centred 
on V'. This KDE step has the effect of smoothing out the variation introduced by using a conditional probabilities to 
estimate a marginal probability, and making the scheme more stable. Increasing the smoothing parameter r results in 
an estimator with a smaller variance, however it also increases the bias of the estimation, so ideally r should be chosen 
to be as small as possible whilst maintaining an acceptable level of variation. 

4.2. Boundary Conditions 

We use the same boundary condition methodology as prescribed by the VRDSMC method, that is for diffusely re¬ 
flecting fully accommodating walls, with temperature T„aii and tangential velocity u„an- Supposing that the Maxwellian 
distribution at the boundary is given by fwaiii'v) = PwaiiP where Pmb is a gaussian probability density, and the 
boundary is the plane x = 0, then the no flux boundary condition is given by 


Pwall 


U J" V:,PMB(y) + j" V^/(V) d\ = 0, 

yr>0 Vv<0 


and similarly for the equilibrium solution 


(37) 


Pwall.e 


, Jv.Pmb, eg(v) dv + J" v^W{\)f{\) d\ - 0. 

yx>0 y;t<0 


(38) 


The second integrals in the above equations are the particle fluxes, and can be estimated by counting the number of 
particles the cross through a wall of area As in a time period At by {\ jAsAf)Nin, and at equilibrium is estimated by 
(1/AsAf) Wi. Also, we can use the analytical properties of the Gaussian distribution to evaluate the first integrals: 


/ 


VxPMBiy) d\ 



(39) 


Therefore, a particle that changes velocity from V to V' when colliding with a wall, changes its weight according to 


yN' = ly(V') = 


/.,(V') 

/(VO 


Pwall,eqPMB,eq(y^ ) 

Pwa iiP mb{'^') 




T^wall 
wall,eq 


Zf" W,- PMB.eg(y') 
Nin Pmb(V0 ’ 


(40) 

(41) 

(42) 


where typically, we choose the temperature of the equilibrium wall boundary condition to be equal to the temperature 
of non-equilibrium wall boundary boundary condition, i.e. T„aU - Twaii.eq- 


5 . Results 

5.7. Homogeneous Relaxation to Equilibrium 

We will demonstrate the effectiveness of this method first with a homogeneous relaxation to equilibrium, i.e. when 
/(f, X, v) = f(t, v) has no spatial component. We start from an initial distribution of particles 


/o(v) = (1/2)(/m(vi;co,co) h-/m(vi;-co,co))/m(v2;0,co)/m(v3;0,co), (43) 
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(a) Without stabilisation (b) With stabilisation 

Figure 1; Homogeneous relaxation towards equilibrium, (a) without KDE, (b) with KDE, smoothing parameter r - 
0.05co. The green squares represent the time-series of the standard MC estimator of the non-equilibrium process; t he 
blue triangles represent the time-series of the standard MC estimator of the process biased to sample from equilibrium; 
the red diamonds represent the time-series of the VREP estimator of non-equilibrium process; the black line represents 
the exact expectation at equilibrium. 


which will relax towards the Maxwellian distribution /m(v, 0 , ^{4/3 )Cq). In figures (lai-(lbi we show how the 
variance reduced estimator performs against the standard estimator, when estimating {|v^using 100 particles, with 
and without the KDE stabilisation procedure. In both cases, the variance of the new estimator is smaller than the 
standard estimator, but the estimator with stabilisation from the KDE reduces the variance even further. 


5.2. Couette Flow 

To test the particle weight variance reduction, we have applied the scheme to sample from a steady-state planar 
Couette flow, and compare to results obtained using a common random number scheme. A Couette flow is a flow 
where the fluid is bounded by two parallel walls moving in opposite directions within their planes, with velocity 
il/waii- For Knudsen numbers Kn - 0.05,0.5,1.0 respectively, Eigures Q, ([^, Q show the variance reduced and 
standard Monte Carlo estimators of the steady-state flow velocity field parallel to the wall, V 2 (-*i), (left) as well as 
the temperature profile across the channel T{x\), for a Couette flow with wall velocity Vwaii = O.Olco, Kn = 0.5, 20 
cells and 100 particles per cell. All the results show a significant improvement of performance over the unweighted 
standard Monte Carlo estimator. 

Next we compare the VREP importance sampling scheme to the CRN correlated equilibrium scheme. Because 
we are interested in the noise of the estimate of the velocity profile across the channel, and the speed of the flow is 
small, we make the simplifying assumption that the steady-state density across the channel is constant. This allows 
us to choose the coefficients A = z/t and D - yllRT^aiilT so that the correlated equilibrium process is distributed 
according to the global Maxwellian, 


/m(x, z, f) = - ^ 3„ exp(—I— ] (44) 

{2nRT„an?l^ \2RT^aiil 

Eigure [^compares the noise-to-signal ratio of the CRN scheme, VREP scheme and standard Monte Carlo estimator 
against signal strength for the samples taken from steady state Couette flow estimator. The results show that for 
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the CRN scheme there is a reduction in the noise-to-signal ratio by a factor of 10 over all tested levels of signal, 
corresponding to a speed up of over 100 times. However, it suffers the same scaling properties with signal size as the 
standard Monte Carlo estimator. In contrast, the importance weighted variance reduced estimator has a noise-to-signal 
ratio that is independent of the signal size as the signal size decreases. This results in an unbounded speed-up over 
the standard Monte-Carlo estimator as the signal size decreases to zero. Because of the independence of the signal 
strength on the noise-to-signal ratio, there is a signal strength where for larger signal strengths, the CRN scheme 
outperforms the particle weight scheme, and for Couette flow we estimate this to be at a Mach number greater than 
0 . 1 . 



xj L xjL 


Figure 2: Couette flow with wall velocity v^aii = O.OIcq, Kn = 0.05, 20 cells and 100 particles per cell. 
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Figure 3; Couette flow with wall velocity Vwaii = O.Olco, Kn = 0.5, 20 cells and 100 particles per cell. 


5.3. Lid-Driven Cavity 

To further demonstrate the effectiveness of the method, we apply it to a lid-driven cavity flow, where the fluid is 
bounded in two dimensions by a square box in the (.r,y) plane, with translational symmetry in the z axis. Three of the 
bounding walls are stationary, and one of the bounding walls moves within its plane at constant velocity t/waii, giving 
rise to a circulatory flow within the cavity. 
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Figure 4: Couette flow with wall velocity Vwaii = O.Olco, Kn = 1.0, 20 cells and 100 particles per cell. 



Figure 5: Comparison of noise-to-signal ratio vs signal size, between standard Monte Carlo, CRN, and our importance 
sampling method (VRFP). 


Figures ([^-(|6b]i show the velocity and non-dimensional temperature field {TITq - 1) of the steady state flow, 
with a lid velocity of t/waii = O.OOIcq for the standard Monte Carlo and variance reduced sampling schemes. The 
results have been averaged over 5000 time-steps, and 10 independent ensembles on a 50 x 50 grid, with an average of 
30 particles per cell. The standard Monte Carlo scheme is not able to pick up the signal, whereas we see clearly that 
the importance sampling scheme is able to recover the signal. In Figure[7]we compute the streamlines of the variance 
reduced flow, and the sheer stress 
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(a) Standard Monte Carlo estimate 



Figure 6: Lid-driven cavity flow. Velocity field and non-dimensional temperature {TITq - 1)/Ma. Kn = 1.0, t/waii = 
O.OOlco, with and without importance sampling variance reduction. 



x/L 


Figure 7: Lid-driven cavity flow. Streamlines and non-dimensional sheer stress 7:\2l{pQRTi)Ms) of the variance re¬ 
duced estimate. 


Figureshows results from lid-driven cavity flows with lid speeds O.lco, O.Olco, O.OOlco and O.OOOIcq. As was 
the case with the Couette flow, the level of noise in each calculation is independent of the lid-speed. 


11 





















0.3 


0 


-0.1 




xxjL 


(a) (/„all - O.lCo 


(b) (/wall = O.OlCo 



Figure 8; Lid-driven cavity flows with different wall-speeds. 50 x 50 grid, 25 particles per cell on average, 5000 
time steps to reach steady-state, thermodynamic fields averaged from 5000 further time steps. The level of noise is 
independent to the wall-speed. 
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6. Conclusion 


In this paper we have developed an importance sampling method for the Fokker-Planck rarefied gas model, that 
assigns weights to each stochastic particle allowing one to sample from an equilibrium distribution. We have demon¬ 
strated its effectiveness in reducing the variance of estimates of thermodynamics quantities for low Mach number 
flows over a range of Knudsen numbers. Significantly, the level of noise in the estimators becomes independent of 
the Mach number for low-speed flows. We believe it to be a versatile and robust method, and because it doesn’t alter 
the basic algorithm of the particle solution scheme, it is able to be used in conjunction with other variance reduction 
schemes such as the CRN method. 
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